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Abstract 

In this paper we present two unconditionally energy stable finite difference schemes for the 
Modified Phase Field Crystal (MPFC) equation, a sixth-order nonlinear damped wave equation, 
of which the purely parabolic Phase Field Crystal (PFC) model can be viewed as a special case. 
The first is a convex splitting scheme based on an appropriate decomposition of the discrete 
energy and is first order accurate in time and second order accurate in space. The second is 
a new, fully second-order scheme that also respects the convex splitting of the energy. Both 
schemes are nonlinear but may be formulated from the gradients of strictly convex, coercive 
functionals. Thus, both are uniquely solvable regardless of the time and space step sizes. The 
schemes are solved by efficient nonlinear multigrid methods. Numerical results are presented 
demonstrating the accuracy, energy stability, efficiency, and practical utility of the schemes. In 
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particular, we show that our muhigrid solvers enjoy optimal, or nearly optimal complexity in 
the solution of the nonlinear schemes. 

Keywords: Nanoscale Materials; Phase Field Crystal; High-order Damped Wave Equation; Con- 
vex Splitting; Finite Difference; Multigrid 



1 Introduction 

Crystalline materials are used in numerous engineering applications. In practice, most crystals have 
nanoscopic imperfections in the form of defects - such as vacancies, grain boundaries, and dislo- 
cations - and controlling, or at least predicting, the formation and evolution of such imperfections 
is a major challenge. Accurate modeling of crystal dynamics, especially defect dynamics, requires 
atomic-scale resolution. The Phase Field Crystal (PFC) equation, a continuum model that has 
attracted significant attention in recent years, has shown promise in this regard O El |T3]. In the 
PFC approach the crystal is described via a continuous phase field (j) that approximates the number 
density of atoms. The field variable (f) admits two qualitatively different solutions, a constant solu- 
tion that represents the liquid phase and a periodic solution that represents the solid phase. The 
local extrema of the periodic solution lie on a near-perfect lattice, mimicking the atomic lattice of 
the crystal. The PFC equation is 

dt(l) = MA{^^ + {l-e)<l) + 2A<l) + A'^^), (1.1) 

where e < 1 and M > 0. This model, which is related to the dynamic density functional theory of 
freezing OE], has the important advantage over many atomistic models that the characteristic 
time is determined by the diffusion time scale and not by that of the atomic vibrations. This 
enables one to capture relatively long time processes, in contrast to atomistic models. The recent 
review by Provatas et al. [13] outlines a wide range of applications of the PFC modeling framework. 

The major disadvantage of the PFC model is that it fails to distinguish between the elastic 
relaxation and diffusion time scales [5l [15] . In [151 [16] , the authors introduced the Modified Phase 
Field Crystal (MPFC) model to overcome this difficulty. The MPFC equation is 

dtt<l) + Pdt(t> = MA + (1 _ + 2A0 + A^^) , (1.2) 



where /3 > 0. The MPFC equation ( 1.2 ) is a nonlinear damped wave equation modeling a viscoelas- 
tic response to perturbations to the density field. In this model, perturbations in the density field 
are transmitted by waves that travel essentially undamped up to a certain length scale determined 
by the parameters. When this length scale is of the order of the size of the system a separation of 
elastic relaxation and diffusion time scales may be practically observed [I5j . 

The MPFC equation, like the PFC equation, is a sixth order evolutive nonlinear partial differen- 
tial equation that cannot be solved analytically in practical circumstances. Efficient and accurate 
methods for approximating the solutions of the MPFC equation are therefore highly desirable. 
However, to date, very little research has been carried out in this direction. 

Because of the close relationship between the MPFC and PFC models, methods for the latter 
equation can be adapted and applied to the former. See, for example, [H |2l [3 ISl [3 121] for some 
recent approximation methods specifically for the PFC model. However, one must take care, as we 
show in the following pages, to adequately account for the wave-like nature of the MPFC solutions 
in the numerical method, especially in the design of provably energy stable schemes. 

Methods specifically designed for the MPFC equation can be found in |10l [161 HBl [E] ■ Ste- 
fanovic et al., [16] employed a semi- implicit finite difference discretization, with a multigrid algo- 
rithm for solving the algebraic equations. They provide no numerical analysis for their scheme. 
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which is significantly different from schemes we propose and analyze. In particular, theirs is not 
expected to be unconditionally solvable or energy stable. The authors did, however, give some 
evidence of the efficiency of their multigrid solver. We will conduct a similar study for our multi- 
grid solver in this paper. The MPFC scheme in [10] is more or less the same as the first-order 
convex-splitting that we devised earlier in |18| [T9] . The first-order convex splitting scheme in our 
work \18\ [19] has two fundamental properties. It is unconditionally energy stable and uncondition- 
ally uniquely solvable. We rigorously analyzed this convex splitting scheme in [THldn], but we did 
not provide a practical solution strategy. That gap is filled here. 

In [51 EI] we presented first and second-order accurate (in time) finite difference schemes for the 
purely parabolic PFC model, based on a convex splitting framework applied to the physical energy. 
The convexity splitting idea - in the context of first order (in time) convex splitting schemes - is 
generally credited to [7]. One of the main advances in [8] - and the more recent work in [13] - was 
the demonstration that the convex splitting framework can be extended to second-order (in time) 
methods as well, in a natural way. While the motivation in [8] was the application to the PFC 
model, the second-order convex splitting idea that we developed is, in fact, rather general. 

The main goal of this paper is to apply the first and second-order convex splitting framework 
to the MPFC equation. The principal challenge in doing so is that the extension to damped wave 
dynamics appears, at first sight, not so straightforward. This obstacle was first surmounted in 
\18\ I19j . as already mentioned, where we extended the first-order convexity splitting idea for the 
MPFC model. Specifically, we showed the existence of a pseudo energy (different from the physical 
energy), which is non-increasing in time for solutions of the MPFC model and to which the usual 
convexity splitting idea can be applied. As a side note, in [18] we used this energy stability approach 
to prove the existence and uniqueness of a global in time smooth solution for the MPFC partial 
differential equation. 

Herein we present, for the first time, a new fully second-order convex splitting scheme for the 
MPFC model. The idea is general and can be applied to other damped wave equations. We also 
compare the first-order convex splitting scheme from [181 with the new second-order convex 
splitting scheme. We propose and test efficient (optimal complexity) nonlinear multigrid solvers 
for both. These two convex splitting schemes are shown to be unconditionally energy stable with 
respect to a pseudo energy, which is different from the physical energy, and mass conserving. Both 
schemes are nonlinear. But both are obtained via gradients of strictly convex, coercive functionals, 
facts which guarantee the unique solvability of the schemes, regardless of the time and space step 
sizes. 

The paper is organized as follows. The MPFC model is recounted in Sec. [2] In Sec. [3] we 
present our first and second-order accurate (in time) convex splitting schemes and analyze their 
basic properties. We present some numerical results in Sec. |4]that give evidence of the convergence 
of the schemes and the efficiency of the multigrid solvers. In that section, we also show a couple of 
practical application of the model. The first is a problem of crystalization, the second, a problem of 
elastic relaxation in a strained crystal. We give some conclusions in Sec. [5} To keep the presentation 
short, we relegate some of the details of the schemes to two Appendices. In Appendix |Aj we give the 
basics of our finite difference discretization of space, and we list some needed summation-by-parts 
formulae. In Appendix [B| we give some details of the nonlinear multigrid solvers that we utilize to 
advance the nonlinear schemes in time. 
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2 The Modified Phase Field Crystal Model 

Consider the energy 

m = l^{\^' + - |V0P + ^(A0)^} rfx, (2.1) 

where : $7 C — t- M is the atomic density field; and a := 1 — e with e < 1. We assume 
il = {0,Lx) X (0, Ly) and (p is J^-periodic. The MPFC equation is the pseudo-gradient flow 

du(l) + f3dt(l) = V • (M(</.)V^) , (2.2) 

where M((f)) > is a mobihty, and 

/i := (50£;(0) = (/>^ + a0 + 2A0 + A^. (2.3) 

The variable fj, is called the chemical potential. For simplicity, we take M to be a constant and 
consider the simpler equation 

dtt^ + ^dtcp = MA/x, fi = cf)^ + a(P + 2A(j) + A'^cj). (2.4) 

Note that in \18\ I19j. we used an alternate version of the equation, namely, 

^du(P + dt^ = Afi, ^ > 0, 

which expresses the model as a perturbed parabolic equation. For the present exposition we will 
abandon this form in favor of Eq. (2.4), since (2.4) is more widely used in the physics literature |15| 

Mi- 
lt turns out that along the solution trajectories of Eq. (2.4) the energy (2.1) may actually in- 



crease in time, on some time intervals. (See Fig.[2]and the corresponding discussion in the Sec. 4.3 



However, solutions of the MPFC equation do dissipate a pseudo energy [l8l[T9|. Also observe that 
Eq. (2.4) is not precisely a mass conservation equation due to the term dtt(f>- However, with the 

introduction of a reasonable initial condition for 5t</>, it is possible to show that / dt(j) dx = for 



all time. 

We shall need a precise definition of the H~J^ inner-product to define an appropriate pseudo 
energy for the MPFC equation. Suppose / G jn G {Q) | /^ndx = 0} =: Lq ($7). Define -0/ £ 
Hp^^ ($7) n Lq ($7) to be the unique solution to the PDF problem 

-A'tpf = f inn. (2.5) 

In this case we write V'/ = — Suppose that f,g £ Lq (f7), then we define 

{f,g)H-^ :=(V0/,V^,)^,. (2.6) 

Note that, via integration by parts, we have the equalities 

U,9)h-^ = - (A-V,5)i2 = - (A-^5,/)i2 = {g,f)H-^ ■ (2-7) 
For every f £ Lq (Q), we define 

= \J{fJ)H-^- (2-8) 
We now recast the MPFC Eq. (2.4) as the following system of equations: 

dtip = MAfi - I3^P , dt(l) =: ip. (2.9) 
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And we introduce the pseudo energy 

f(<A,V):=i^(0) + ^IIV'llH-M (2-10) 

which requires that J^ip <hi. = 0. This is the case as long as the initial condition dt(j){x, 0) dx = 
is satisfied \18i. il9j . In what fohows, we will use the initial data 

^(x,0) = 9t0(x,O) EE in Q, (2.11) 

so that dt(l){x, 0) dx = is trivially satisfied. In any case, as long as ip = dt(p is of mean zero, a 
formal calculation shows that (see ^51 |T9] ) 

dt£ = iii, dt4>)^, + ^ (V^, 9t^)^-i = ^ 0' (2.12) 

which guarantees that the pseudo energy is non-increasing in time. 

At this point we observe that the energy S((/>) admits a splitting into purely convex and concave 
energies, given hy E = Ec — Eg, where Ec and E^ are both convex functionals. In particular, we 
have the convex splitting of the form 

Ec = \ M% + f ||0|li. + I \m\l, , E, = \\Vm, . (2.13) 

This observation is the key toward the construction of a (first-order) convex splitting scheme, as 
was demonstrated in \18\ [T9] . It will also be the basis of the new fully second-order scheme that 
we present in the following sections. 



3 Finite Difference Schemes and Their Properties 



In this section we present two finite difference methods for the MPFC equations. The first is a 
first-order accurate convex splitting scheme that was introduced and analyzed in |18t I19j . The 
second is a new, fully second-order accurate convex splitting scheme for the MPFC equation that is 
similar to a scheme we presented earlier for the PFC equation [S] . In the following we will show that 
both schemes are discretely mass conserving and unconditionally energy stable. In addition, both 
schemes are unconditionally uniquely solvable. We refer the reader to App.[A]and also [20^,21] for 
a description of the finite difference framework that we employ to define and analyze our schemes. 

We begin by defining discrete versions of the energy (2.1) and pseudo energy (2.10). For the 
finite difference grid function (f> £ Cmxn define 

1 



F{<P) :-- 



4 a 
4+ 2 



V 



2 1 

2 + 



For 



4 " 2 " " 2 

£ Cffixn and ip G H, define the discrete pseudo energy as 

1 



l^ 



(3.1) 



2M 



-1 ' 



-1' 



(3.2) 

norm is a finite-dimensional 



where H is the space of mean zero cell centered functions and the 
equivalent to the continuous norm. See App. \a\ in particular (A. 14) - (A. 18), and also [19J. 
Note that the discrete energy F admits a convex splitting. In particular, if the grid function 
</> £ Cffixn is periodic, then the energies 

^ ^ " ^ " - ' (3.3) 



F, 



4 a 
4+ 2 



2 1 
2+0 



\l 



F,, 



are convex 

-|- a 



and F = Fc — Fa. Furthermore, the discrete variations (gradients) are 6(j,Fc{<j)) 



^^-'^ + Alcl) and 6^F,{(t)) 
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3.1 First-Order Convex Splitting Scheme 



The first-order convex-splitting scheme, introduced in [THl [T9|, is constructed by respecting the 
convexity splitting F = — F^. and is defined as foUows: find (j)'''^^ , ■ip'''^^ , fi'''^^ , k'^'^^ G Cmxn 
periodic such that 



kk+l 



5rhF,, 



k+l 



(3.4) 
(3.5) 

(3.6) 
(3.7) 



where s > is the time step size, which we wiU assume is unchanging, and ip^ = 0. Notice that the 
part coming from the convex energy, Fc, is treated implicitly, and the part coming from F^, is treated 
explicitly. In this case, such a treatment leads to a nonlinear scheme, though the nonlinearity is 



quite mild. This scheme can be simplified by using (3.5) to eliminate the unknown ij)^'^'^ from (3.4) 
Therefore, an equivalent form is 



sM 



k+l 



Ak+1 



k+l 



+ 



/3 + 



(3.8) 
(3.9) 



Note that Eqs. (3.8) and (3.9) decouple, and ip^'^^ can be obtained after cj)^^^ is known. The three 



equations (3.6), (3.7), and (3.8) are solved simultaneously using the non-linear multigrid method 



that is described in App. B.l Then tp'''^^ is updated via (3.9). 



3.2 Second-Order Convex Splitting Scheme 

Following our work in [U [H] we propose the following second-order convex splitting scheme: given 
c/)'^, i/j'^, G Cmxn periodic, find 4>^~^^ , ip^^'' , g Cfnxn periodic, such that 



- ^p^ 



X 

Ah 



Lk+\ 



(3.10) 
(3.11) 
(3.12) 

(3.13) 



where s > is the time step - which, again, we take to be unchanging for simplicity - and 



b'^ + ip^ 

2 ' 



^fc+i _^ ^k 
2 ' 



kk+l 



+ 



Ik+i 



3/' - (p^-^ 
2 



We enforce the initial conditions ip^ = {) and (p^^ = (p^ . Note that the approximation (p^^ = cp^ will 
result in a first order (in time) local truncation error at the first step of the scheme. This does not, 
however spoil the global second-order convergence of the method, as our tests later show. 

Here again, we respect the convexity splitting noted earlier. The convex contribution to the 
chemical potential is treated implicitly, but now using a second-order secant type approach as 
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in [HlITlj. The concave contribution is treated explicitly using a second-order extrapolation. This 
leads to a nonlinear, two-step method. As before, we can decouple some of the equations in the 
scheme. In particular, eliminating ip'''^^ from (3.10) we obtain 



fc+1 



sM 

2 

s 



tfc+l 



(3.14) 
(3.15) 



so that ^^^^ may be calculated after (/)^"^^ is known. The three equations (3.12), (3.13), and (3.14) 



are solved simultaneously using the non-linear multigrid method that is described in App. B.2[ 
Then ■0'^"'"^ is updated via (3.15). 



3.3 Mass Conservation and Unconditional Unique Solvability 

The above schemes are discretely mass conserving and uniquely solvable for any time step s > 0. 



These facts were proven for the first-order convex splitting scheme (3.6) - (3.9) in [19j. We now 



provide the proofs for the second-order convex splitting scheme (3.12) - (3.15) 



Theorem 3.1. (Mass Conservation): The first- order convex splitting scheme (3.6) - (3.9) and the 
second-order convex- splitting scheme \3.12 ) - (3.15) for the MPFC equations are mass conserving 
for any time step s > provided solutions exists to the respective schemes. 

Proof. The details for the first-order scheme are contained in IIQ, Thm . 4. 2]. S uppose that 



'f/;'^"'"^, /i'^"'"2 J K*^+2 J is a solution to the second order scheme (3.12) - (3.15). Summing 



Eq. (3.10) and using summation-by-parts, specifically Prop. A. 2, we have 



fc+1 



sM Ahfi' 



sp 



-sM 



-s/3 



^fc+l ^ 



sM 



Dyl 



sf3 



(3.16) 



This gives the relation 



fc+1 



2 



1 + 



13s 



This fact, together with the initial condition i()^ = 0, ensures that (V'^Hl) 



(3.17) 

for all A; > 0. Now, 



observe that from (3.11) 



and the result follows: (fl!)'^"'"^||ll 



if and only if 



^fc+l _^ 



0, 



(3.18) 



□ 



At this point it is worth emphasizing that for the simplified initial condition -i/;'^ = we have 
= for all > 0. In other words, is a mean zero, periodic function for all /c > 0. It 
would be enough for our purposes to enforce a more general initial condition, namely, (V'^||l) = 0, 
though we do not pursue this level of generality here. 
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Theorem 3.2. (Unconditional Unique Solvability): The first-order convex splitting scheme (3.6) 



(3.9) and the second-order convex splitting scheme (3.12) - (3.15) for the MPFC equations are 



uniquely solvable for any time step size s > 0. 

Proof. The unconditional unique solvabihty for the first-order convex splitting scheme was proved 
in [i9i Thm. 4.3], and we skip the details for that case here. The proof for the second-order convex 



splitting scheme (3.12) - (3.15) closely follows our techniques from [19] for the first-order MPFC 



scheme and [8] for the second-order PFC scheme. Specifically, we first construct the functional 



where 



and 



1 

mn 



G{4>) ■■= , 

b^\\l), 



2 

-1 



+ 



sM 



R{(t>), 



+ 2A. 



(3.19) 

(3.20) 
(3.21) 

(3.22) 

This, together with the properties of the "—1" inner product and norm given in App. [Aj implies 
that G{(j)) is strictly convex and coercive over the set of admissible functions, the hyperplane 



Q(0) : = 



4 



+ 



1 + 



a 



Following arguments in [8], we see that -R(</>) is strictly convex and 



A={c 



1^ and <j) is periodicj , 



and its unique minimizer 



kk+l 



€ A satisfies the discrete Euler-Lagrange equation 



5^G ( c/^'^+i ) = -A 



-1 
h 



/3 + 



+ 



sM 



/+! + c = 0, 



where C is a constant. This is equivalent to 



sM 



/3 + 



(3.23) 



(3.24) 



which is just Eq. (3.14). Thus minimizing the strictly convex functional (?((/>) over the set of 



admissible functions A is the same as solving the second-order convex splitting scheme (3.12) 



(3.15). This completes the proof. 



□ 



3.4 Unconditional Energy Stability 

The following theorem was proven in |19] . 

Theorem 3.3. (Unconditional Energy Stability of the First-Order Scheme): The first- order convex 



splitting scheme (3.6) - {3.i^ is unconditionally (strongly) pseudo energy stable with respect to (3.2) 



meaning that for any time step s > 0, 



(3.25) 
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We now proceed to show a similar result for the new second-order convex splitting scheme. 

Theorem 3.4. (Unconditional Energy Stability of the Second-Order Convex Splitting Scheme): The 
second-order convex splitting scheme (3.12) ~ (3.15) is unconditionally (strongly) energy stable with 
respect to the discrete energy 



In other words, for any s > and any h > 0, 



kk j.k 



-fc-1 



More specifically, 



+ s 



M 



where 



Proof. Using the identity 



2 Uk+1 



D 



2A, 



2 



l2 i/fc 



+ 



2 1 

. + 2 

2 1 



2 2 



l,k+l _ A.k 

2 
2 



and Eq. (3.12) we obtain 



F 



k+i _ 

1 



+ 



-F 
1 



2 
1 

2 

V^iAk+i 



kk+l 



kk-l 



2<P' + 



Next, using Eq. (3.10) and properties of the norm || • ||_^ (see App.|Aj), we have 



1 

2M 



k+l 



2 1 

-1 ~ 2M 



M 
M 



sAhfi' 



k+k 



k+k 



-1 



s/3 






s/3 







Adding Eq. (3.31) and Eq. (3.32) we have the result. 



(3.26) 
(3.27) 

(3.28) 
(3.29) 



(3.30) 



(3.31) 



(3.32) 

□ 
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4 Numerical Results 



In this section we present some numerical results for the convex splitting schemes that illustrate 
the convergence rates and energy stability of the schemes. We conclude the section by applying 
the model to a problem of crystallization in an undercooled melt. 



4.1 Convergence Tests 

In order to study the theoretical convergence of the two schemes we evolve the same initial data 
in time with increasingly finer grid spacing and compare the solutions. We choose a domain 
n = (0, 32) X (0, 32) with an initial density field 



ix,y,t = 0) 



0.07-0.02 cos 



27r(x - 12) 
32 



sm 



( My 



1] 



V 32 



2/'vr(x + 10)\ 2/^7r(y + 3) 
+0.02 cos^ I -^^t:; 1 cos^ ' ^ 



32 



32 



-0.01 sin^ 



Attx 



. 2 /47r(y-6) 

, sm 

32 y V 32 



(4.1) 



The parameters are taken to be M = 1, e = 0.025 (a = 1 — e) and /3 = 0.9, and we evolve the 
system to final time tf = 10. For these initial data and parameters, the system evolves toward a 
nonuniform steady state. 

The time step refinement paths are given by s = 0.025/i^ for the first-order convex splitting 



scheme (3.6) - (3.9) and s = O.Obh for the second-order convex splitting scheme (3.12) - (3.15) 



Note that these refinement paths have nothing to do with CFL-type stability constraints. The 
schemes are unconditionally stable. The important point is that, using these paths, in both cases 
the global error is predicted to be £{tf) = 0{h'^), and this will be verified in the tests. The fine 
details of the convergence tests are the same as those described in [HI [U] . We skip the specifics 
and direct the interested reader to those references for a more complete discussion. The essential 
idea is to compare solutions at successively finer resolutions, one representing the coarse resolution, 
h = he, and one the fine resolution, h = h 



and each one half the size of 



QQ QQ oQ on 

The mesh spacings are taken to be /i = gj) ils' ^ 
the previous. The results of the tests are given in Tabs. [T] and [2] for the first and second order 
convex splitting schemes, respectively. The test data provide evidence that both schemes are 



second order accurate in space. They suggest that the scheme (3.6) - (3.9) is first-order accurate 
in time, i.e., £i{tf) = 0{s) + 0{h?'), while (3.12) - (3.15) is second order accurate in time, i.e., 
S^itf) = 0{s^) + 0{h'^), as expected. 



4.2 Multigrid Efficiency Test 

We now present a test, similar to one from [20j . that shows the efficiency of our nonlinear multigrid 
solvers described in App. [Bj Specifically, we fix the time step size at the relatively large value of 
s = 10.0 and vary the spatial step sizes: ^ = ||) g|) and We use the initial data 

(4.1), and the parameters are the same as in the last test: M = 1, e = 0.025 (a = 1 — e), /? = 0.9, 
Cl = (0,32) X (0,32). We advance the solution using our second-order convex splitting scheme 
(3.12) - (3.15) by exactly 20 time steps, so that the final time is tf = 200. We use -^max = 1 in the 
definition of the smoothing operator (cf. App. B.2). In Fig. [T]we report the reduction of the norm 
of the residual for each v-cycle iteration and for each space step size h. The reader will observe that, 
independent of h, the norm of the residual is decreased by a factor of about 5.5. This indicates that 
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the norm of the error - not the discretization error, but the algebraic iteration error of our imphcit 
scheme - is reduced by roughly the same amount. This gives strong evidence that our multigrid 
solvers have optimal, or nearly optimal, complexity. The results of the corresponding test for the 



first-order scheme (3.6) - (3.9) are expected to be similar; we skip this test for the sake of brevity. 



4.3 Energy Stability Test 

In this section, we demonstrate the energy stability of our second-order convex-splitting scheme 
( |3.12 ) - (3.15). We again use the initial data (4.1). The physical parameters are M = 1, e = 0.025 
(a = 1 — e), = (0,32) x (0,32), and /3 = 0.01. These are the same as in the last couple of 
tests, except for the smaller value of /?, which changes the time scaling and also emphasizes the 
wave nature of the equation more than in the other tests herein. We use h = ^ and s = 0.1 and 
advance the numerical solution to tf = 100. We plot the calculated discrete energies J-, defined in 
Eq. (3.2); defined in Eq. (3.26); and F, defined in Eq. (3.1), against time in Fig. 2 As pointed 



out earlier, the energy F may be increasing in time; this phenomenon is observed here. This test 
is somewhat contrived, since the value of f3 may be smaller than what is typically used in practice. 
If (3 is large, in fact F is often observed to be non-increasing. However, this shows the point that 
F is the incorrect energy for this dissipative system. The only energy that is guaranteed to be 
non- increasing is J-"; this property is observed here. But, note that J- is nearly identical to This 
phenomenon is expected when s is sufficiently small. 



4.4 Analysis of Effective Time Step Sizes 

We now proceed to study the effect of the choice of time step s on the predicted crystal dynamics 
following similar tests in [3l [8]. We perform simulations on a domain of size = (0, 128) x (0, 128) 
and initial data (pij = (j) + ijij where (p = 0.07 and r]ij is a uniformly distributed random number 
satisfying \rjij\ < 0.07. The parameters are M = 1, (3 = 0.9, e = 0.025 (q = 1 — e) and h = 1. 



To create a baseline, we first perform a simulation with the fully second order scheme (3.12) - 
(3.15) with the relatively small time step size s = si = 0.01. The density fields at time t = 350, t = 
500, t = 1500 and t = 2500 are shown in the first row of Fig. [3} Then we perform successive 
simulations with S2 = 1, S3 = 10 and S4 = 20 (second, third, and fourth columns of Fig. [3]) again 



using the second-order convex splitting scheme (3.12 ) - (3.15 ), with all other physical and numerical 



parameters being the same. Following the procedure in [U [8] we compare frames of the different 



simulations when the pseudo energy T (3.2) values are the same, rather than when the times are 
the same. Each column of Fig. [3] represents configurations of (f> that have the same pseudo energy 
T. The scaled pseudo energy is plotted for the s = si = 0.01 case in Fig. |4j it is observed to be 
non-increasing in time. 

(Islel) 



Now, we repeat the test using the first-order convex splitting scheme ( |3.6[ ) - (3.9) and report 
the results in Fig. [5j But here, for the first row of Fig. [5| we plot (j) at the times for which the 
pseudo energy values match those for (p represented in the first row of Fig. [3j calculated using the 
more accurate second-order scheme (3.12) - (3.15). Again, the <j) configurations in each column of 
Fig.[5]have the same pseudo energy values, which explains the discrepancies in the corresponding 
times. 

Inspecting Figs. [3] and [Sj the qualitative agreement of the density fields (j) corresponding to 
the same pseudo energy values (plots in the same respective columns) is apparent. However, 
the numerical times required to reach the specific pseudo energy levels can be wildly different. For 
example, for the first-order convex splitting scheme for the S4 case (Fig.jsjrow 4) it takes t = 122640 
to reach the energy corresponding to the second order scheme at t = 2500 for si, about 50 times 
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longer. Meanwhile, the second order scheme reaches the same energy level for 54 case at t = 5520 
indicating the higher time accuracy of the second order scheme. The second-order convex splitting 
scheme is clearly more accurate, at least at a qualitative level, for this practical test. 

Following [31 [8], in order to study the quantitative difference in respective solutions (obtained 
with different time step sizes s) at the same final time tf, we define the scaled difference 

^ _ ll'/'(-,t/;si) -(/>(-, t/;s) 11^ 

Ui-,tr,s,)\\, ' ^'-'^ 

where we choose tf = 350. The other parameters are the same as in the previous test: M = 1, 
/? = 0.9, e = 0.025 (a = 1 — e) and h = 1. Figure [6] shows the log-log plot of the scaled difference 
obtained using the second order convex splitting scheme at time t/ = 350 and the baseline time 
step size s = si = 0.01. The results are quite similar to that observed for the PFC schemes in [3l[8j. 
The accuracy is observed to increase as the time step is decreased, while the difference (the error) 
tends to saturate for large time steps. 

4.5 Growth of a Polycrystal from a Supercooled Liquid 

In this test, we apply our second-order scheme to solve a problem of polycrystallization in a su- 
percooled liquid. We start our simulation with a constant density field (j) = 0.285 on a domain 
Q = (0,804) X (0,804). We place three random perturbations (small in spatial extent) along 
the bottom y = line of the domain as seeds for nucleation. The domain is discretized using a 
2048 X 2048 grid, i.e., h = 804/2048.The parameters are set to be M = 1, e = 0.25 (a = 1 - e), 



and /3 = 0.9, and snapshots of the crystal microstructure are shown in Figs. 7a - |7d[ 

Here the boundary conditions are periodic in the horizontal direction and homogeneous Neu- 
mann on the top and bottom boundaries. Our theory can be easily extended for these conditions. 
This case study serves to show the applicability of the model to a physical problem; it also shows 
the relative flexibility of our multigrid/finite-difference framework over spectral and pseudo-spectral 
methods. In particular, it is a triviality to modify the nonlinear multigrid solver to adjust to the 
present mixed boundary conditions from fully periodic boundary conditions. This is not the case 
for spectral methods, where the entire method mat have to be redesigned for different boundary 
conditions. 

Three different crystal grains are observed to nucleate from the seed sites and grow. These 
grains eventually grow large enough to form grain boundaries as they impinge upon one another. 



In Figs. 7b - 7d, the grain in the middle forms boundaries with the grains on the left and right. 
However, the orientations of the two crystallites on the right and left extremes are so similar to one 
another that they relax without a definitive grain boundary, but form point defects (about x = 50, 
t = 1000 and t = 2000) along where the grains impinge. 

In order to study the effect of the damping term we repeat the simulation with the same initial 
data and parameters, with the exception that we take /3 = 10 to obtain higher relative damping 
than in the last test. A single snapshot at time t = 3000 for this new simulation is shown in Fig.|8j 
The increased damping manifests itself in two significant ways compared to the previous case. The 
first is in the qualitative features of the grain boundaries, and the second is in the crystal growth 
rate. 

Some of the second can be explained by the change in the diffusion time scale with the modi- 
fication of /3. If we rescale the present case so that the diffusion time scale is the same as that of 
the last, we obtain 

^] dss<p + 0.9ds<P = Afx, (4.3) 
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where s = j^t. Thus t = 3000 corresponds to the synchronized diffusion time s = 270. The earhest 



snapshot shown from the last simulation is at time t = 500 (Fig. 7a). Clearly for the high damping 
case the crystal is growing into its melt faster (from the perspective of the synchronized diffusion 
time scale) than in the low damping case. 

On the first point, the increased damping (/3 = 10) of the propagation of waves allows the 
left grain and the middle grain to merge into one crystallite of a single grain with a single point 
defect located at the bottom (at x = 250) unlike the case with /3 = 0.9. One should compare the 
microstructure with that of Fig. |7b| {t = 1000). More quantitative differences between the low and 
high damping cases can be observed from the relative pseudo energies, which are plotted for both 
simulations in Fig. [9j 

4.6 Effect of Applied Strain on Coherent Crystal 

The key difference between the the PFC and MPFC models is the ability of the latter to capture 
the difference between the "elastic" relaxation and "diffusion" time scales [15]. To show the power 
of the MPFC model we perform a test that highlights this difference. Specifically, we will study the 
relaxation of a strained crystal by tracking the displacement field with respect to the strain-free 
crystal configuration over time. 

We take e = 0.6 (a = 1 — e), M = 15^, h = j^^, s = 0.025/i^, and use the second-order convex 
splitting scheme for the tests. We use two different values of (3 in the tests. The grid sizes are 
m = 512 (grid points in x direction) and n = 400 (grid points in y direction). A solid periodic 



crystal, 20 atomic layers thick, was formed in a liquid bath, as shown in Fig. 10 in the coexistence 
region, i.e., where the solid and liquid coexist at equilibrium. The coexistence solid and liquid 
densities for e = 0.6 are (j)s = 0.395 and 4>i = 0.57, respectively [5]. The solid region density field is 
chosen to be the one given by the single mode approximation to the equilibrium density [5]: 



4>{x,y) = (t)s + A 

where 



cos (qtx) cos —p= cos — ^ 



(4.4) 



A=Us + ^Vl5^^^s, Qt = ^- (4.5) 
5 15 2 

The bottom row of atoms in this crystal were pinned in this configuration and then annealed 

to equilibrium, leading to the strain-free microstructure shown in Fig. [9j Then a shear strain was 

applied by pinning the density field of the atoms along the top row and shifting the density field 

of the bottom row of atoms by 4% of the horizontal lattice spacing, i.e., 4% of 27r, to the left. 

The pinning of atoms and the application of strain is accomplished by following the technique of 

Stefanovic et al. [151 [16]. The idea is to modify the continuous free energy, E((j)), to obtain 

E{<l)) := E{<j)) + J^M{x,y) (^(f>{x,y) - 4>ix,y)y dxdy, (4.6) 

where > is a weight function and is a prescribed crystal configuration. We take = 2 in 
horizontal strips of vertical thickness vr (about half the lattice spacing) around the top and bottom 



of the crystal, and otherwise A4 is set to zero. The function (j) is chosen similar to (4.4) but modified 
slightly so as to effect the desired pinning and shear [TB] . The incorporation of the modified 
energy into our second-order convex splitting scheme is straightforward. We simply replace ^^'^'2 
with 

/i'^+l +2 ((^i.'^+i -(^)X. (4.7) 
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The horizontal strain at an atom position is defined VI cl Sxx 

:= -FT^, where Ax is the horizontal 
displacement measured from the reference configuration. At steady state, we expect the strain to 
vary linearly from 4% to 0% along the height of the crystal. Figure 11 shows the time evolution 
of the strain field for two computational cases: f3 = 0.2 (triangles in Fig. 11, low damping case) 
and (3 = 9 (squares in Fig. 11 high damping case). Note the oscillatory nature of the low damping 
MPFC case (triangles). In this case the crystal rapidly relaxes to the strained equilibrium state, 
represented by the straight line. The high damping MPFC case behaves much more like the PFC 
model, as expected, with a slow relaxation to equilibrium, but without the oscillatory overshoot 
(even at long times, not shown) that is seen in the low damping case. 



5 Concluding Remarks 

In this paper we have presented two unconditionally (pseudo) energy stable finite difference methods 
for the MPFC model. Both are based on convex splittings of the physical energy. One is first-order 
accurate in time and appeared in our previous works [181 HH] • The other is second-order accurate 
in time; and the way we deal with the the chemical potential term in the scheme is similar to that 
in our second-order convex splitting scheme for the PFC model [8j. What makes this a non-trivial 
extension of the previous method is the treatment of the second-order time derivative dtt4>- The 
second-order approximation of this term relies on the introduction of a new variable ip := dt(j). 

We provided numerical evidence of the respective orders of accuracy and convergences of the two 
schemes, the efficiency of our nonlinear multigrid solvers, and the (pseudo) energy stability of our 
schemes. In particular, we showed that, even though our schemes are nonlinear, they can be solved 
with optimal (or near optimal) complexity using nonlinear FAS multigrid solvers. We demonstrated 
that the second order accurate scheme generally captures the time evolution with greater accuracy 
and efficiency than the first order scheme, which is of course what would be expected. We concluded 
the paper with a couple of practical computational examples, one of solidification and one dealing 
with elastic relaxation in a strained crystal. 

We take this opportunity to emphasize the significance of the second-order convex splitting 
scheme and attempt to put the work in some broader context. The novelty of our second-order 
convex splitting scheme is not that it enjoys an unconditional pseudo energy stability (ES) property 
and not that it enjoys an unconditionally unique solvability (US) property but that it enjoys both 
properties simultaneously. One can always construct a scheme that satisfies one or other of the 
two individually. For example, fully explicit schemes are usually always unconditionally uniquely 
solvable (trivially). And one can always construct a second-order energy stable scheme using a 
secant-type (or Crank-Nicolson-like) approach. Here, roughly speaking, given an energy E, one 
takes 

._ E {4>'^') - E (0^) 

^ - ^k+i _ ' y^-^) 

as in P]. (In fact, this is the idea we use to deal with the convex portion of the energy (and the 
chemical potential) in our own second-order schemes, here and in [^.) While, this will lead to 
an unconditionally energy stable scheme, the question of solvability is much harder. Typically, 
unconditional solvability must be sacrificed. 

As we have pointed out, our scheme enjoys both properties ES and US. But this comes at a 
price. Namely, we use a multistep method, and it is more difficult in this setting to employ adaptive 
time stepping. What would be ideal is a one-step, second-order scheme enjoying both properties 
ES and US. In this case time adaptivity could be more straightforward. Regardless of whether it is 
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using a new second-order scheme or the present one, adaptive time and space stepping strategies 
constitute an important future research direction. 
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A Finite Difference Discretization on a Staggered Grid 

For simplicity, we assume that Q = (0, L^) x (0,Ly). The framework that we describe has a 
straightforward extension to three space dimensions. Here we use the notation and results for 
cell-centered functions from [20\ [2T] , see also [8l [19] . The reader is directed to those references 
for more complete details. We begin with definitions of grid functions and difference operators 
needed for our discretization of two-dimensional space. Let = (0, Lx) x (0, Ly), with Lx = m ■ h 
and Ly = n ■ h, where m and n are positive integers and /i > is the spatial step size. Define 
Pr := {r — ^) ■ h, where r takes on integer and half-integer values. For any positive integer £, define 
Eg = {pr \ r = ^, . . . ,£ + 1} , Ci = {pr \ r = 1, . . .,£}, Cj = {pr ■ h \ r = 0, ...,£+ 1}. Define the 
function spaces 

Cmxn = W-.CmXCn^R}, C^xn = X ^ M} , (A.l) 

Cmxn = {0 : Cm X Cn — ?• M} , Cmxn = {0 : Cm X Cfi — )• M} , (A. 2) 

£mxn = {u:EmXCn^m, C'xn = {v. Cm X En ^ , (A.3) 
£^^^ = {u:EmXCn^m, £mxn = {v. Cm X E^ ^ R} . (A.4) 

We use the notation (pij := <p{pi,Pj) for cell-centered functions, those in the spaces Cmxn^ ^mxrit 
Cmxn, or Cmxn- In Component form east-west edge- centered functions, those in the spaces £mxn 
or £mxn-: identified via ^^j+i j := u{p-_^_i,pj). In component form north-south edge-centered 
functions, those in the spaces £mxni ^mxni ^"^^ identified via u^_^i ■ := u{p-^i,pj). The functions 
of Vmxn are called vertex- centered functions. 

We need the weighted 2D grid inner-products ( • |l " )' [ " II " lew I " II ' Ins ^^^^ are defined in [201 [2T]. 
In addition to these, we also need the following one-dimensional inner-products: 

m n 



2 u ' I 2 ' 



where the first is defined for f,g^ ^mxni and the second for f,g^ ^mxn- 

The reader is referred to [20| [2T] for the precise definitions of the edge-to-center difference 
operators dx : £mxn ~^ Cmxn and dy : £mxn ~^ Cmxn] the X— dimension center-to-edge average and 
difference operators, respectively. Ax-, Dx : Cmxn — ^ ^mxn'^ the y— dimension center-to-edge average 
and difference operators, respectively, Ay, Dy : Cmxn ^mxn'i the standard 2D discrete 
Laplacian, A/j : Cmxn ~^ Cmxn- 

In this paper we use grid functions satisfying either periodic or homogeneous Neumann boundary 
conditions, or some appropriate combination of these. Specifically, we shall say the cell-centered 
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function (p € Cfnxn satisfies homogeneous Neumann boundary conditions if and only if 



4>i,n+l 



4'i,ni 



hfl = (f>i,i, i = 0, . . . ,m + 1 . 



(A.6) 
(A.7) 



We say that the cell-centered function (p G Cmxn satisfies periodic boundary conditions if and only 
if 



4'i,n+l 



Pifi 



'^m,j ) j — 1, . . . , n , 
bi„, i = 0, . . . ,m + 1 . 



(A.8) 
(A.9) 



We will use the grid function norms defined in [2^ [21] . The reader is referred to those works 
for the precise definitions of || • Hg, || • \\^, \\ • ||p (1 < p < oo), || • ||q 2, || • ||i 2; and \\(j)\\2 2- 

Using the definitions given in this Appendix and in [2^ ^IT\ , we obtain the following summation- 
by-parts formulae: 

Proposition A.l. (summation-by-parts) If (j) £ Cmxn U Cmxn and f G S^xn then 



and if (f) e Cmxn U Cfnxn and f e£, 



-h (A 



14./) 



fm.4- 4 . 



mxn then 



-h^ mdyf) 

-h (aU 



Proposition A. 2. (discrete Green's first identity) Let (p, ip £ Cmxn- Then 



|A^V) 



h [A^di 



+ h A^ 



h [Ay(P^l DytP^l^ +h {Ay(P^^^^ Dy^^ 



Proposition A. 3. (discrete Green's second identity) Let (p, ip G Cmxn- Then 



h^ {Ah 
+h (a 



Axipi 



+h (A 



Dyt 



y'fi.,n+l 



h D 



'yH'^^n+i 



Ayi^_ 



y'Ti.^n+i 



-h [A 



(A.IO) 



(A.ll) 



(A.12) 



(A.13) 



All of the boundary sums above vanish when the boundary conditions are periodic or homoge- 
neous Neumann. 
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We need the definition of the following space in the analysis of our scheme: 

H:= {</.GC„xn|(</'||l) =0}. 
We equip this space with the bilinear form 

for any (pi, (j)2 & H, where ^pi S Cfnxn is the unique solution to 

- Ahipi = (pi, ipi periodic, (■0i||l) = 0. 
The proof of the following can be found in [19j. 

Proposition A. 4. ((/>i||</>2)_x is an inner product on the space H. Moreover, 



h 'Pill'Ps 



)• 



Thus 



defines a norm on H . 



(A.14) 



(A.15) 



(A.16) 



(A.17) 
(A.18) 



B Nonlinear Multigrid Solvers 

In this Appendix we provide the details of the nonlinear multigrid algorithms that we use for time 
stepping the semi-implicit numerical schemes proposed in previous sections. We note that the 
first-order convex splitting scheme was first presented in |19], but no implementation strategy or 
simulations were presented. 

B.l Multigrid Method for the First-Order Convex SpHtting Scheme 

For the first-order convex splitting scheme (3.6) - (3.9), given (p^ , ijj^ ^ Cmxn, we need to find the 
grid functions (p, ip, ii, k G Cmxn satisfying 



0- 



( sM \ 



/3 + 



-lb''- 

1 'f^^J' 



(p% + a(pij + 2A/,<; 
^h<Pij, 



(B.l) 

(B.2) 
(B.3) 

(B.4) 



where we have dropped the superscript k + 1. We solve this system by observing that the first 



three equations decouple from the fourth. Thus, we solve for (p, fi, k £ Cmxn using (B.l) - (B.3) 
and then update tp G Cmxn using (B.4). 



Equations (B.l) - (B.3) are nonlinear and can be solved using a nonlinear multigrid algorithm 

The first task it to identify a nonlinear operator 



similar to the one proposed in [8j. See also 
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N and a source term S, such that N = S is equivalent to (B.l) - (B.3). Let u = {(f>, jj, k)'^ . The 
3 X m X n nonhnear operator N(u) = (A^(-'-)(u), A^'^^)(u), A^('^^(u))-^ is defined via 



sM 



f{u) := ^,,-4 
N^f\u) := Kij-Ah^,j. (B.7) 



^ij i'^) ■= fJ-ij - <Pij - (^(Pij - ^hKij, (B.6) 

r(3)/ 



The 3 X m X n source S= is defined via 



Clearly the system ( |B.l ) - (B.3) is equivalent to N(u) = S 



The system N(u) = S can be efficiently solved using a nonlinear FAS multigrid method. Ex- 
tensive literature is available and we refer the reader to |17j for a detailed review. Since we are 
using a standard FAS V-cycle approach, as in [U |9l [20] we only provide the details concerning how 
the nonlinear smoothing operation is performed. The other components of the algorithm can be 
constructed by the interested reader by consulting [17J. 

For the smoothing operator we use a nonlinear red-black Gauss-Siedel method. Let i denote 
the smoothing iteration. By ^max we mean the maximum number of smoothing iterations. We 
point out that the nonlinearity in the operator N comes only from the cubic term . This can be 
handled with a local linearization. We commonly use either a local Newton-type linearization or a 
simpler local Picard-type linearization, respectively, 

3(0!,)' or (ef (m 

To simplify the description of the smoother, we write the details assuming a lexicographic ordering is 
followed and we employ the Picard linearization. To get a proper description of the implementation 
of the red-black ordered Gauss-Seidel method, which is what we use in our codes, we refer the reader 
to |17j . The smoothing scheme can thus be written as 

(B.IO) 

^^^' + 4^' = sS + ^U,^, + 4^^, + ^l + <Pl^.y (B.12) 



At each grid cell we solve for the unknowns ^^j^^, using Cramer's rule. One 

complete smoothing pass is effected when all of the grid points have been visited exactly once. 
The smoothing operation is complete when ^max full smoothing passes (or smoothing iterations) 
have been conducted. 



18 



B.2 Multigrid Method for the Second-Order Convex SpHtting Scheme 

In the case of the fuhy second order scheme (3.12) - (3.15), given (p^ ^ijj^ ,(f)^^^ G Cmxn we need to 
find the periodic grid functions ip, fi, k £ Cmxn satisfying 

sM . 2 



/3 + 



(B.13) 



(B.15) 
(B.16) 



where we have dropped the superscripts on the unknowns for the sake of brevity. As before, we 
solve for ^, k G Cmxn using (B.13) - (B.15) first and then update ip G Cmxn using (B.16). 

With the same notation as before, the nonhnear operator N(u) = (A^(^''(u), A^(^^(u), A^^'^)(u))"^ 
is defined component-wise via 







:= 4>ii 


Ms 

^ _j_ 2 '-^hfJ'tj J 








'■= ^J'ij 










'■= Kij 







(B.17) 

(B.18) 
(B.19) 



and the source term S is defined component- wise via 



S. 



(1) ._ 



+ 



(2) 



^(3) ._ 1 A rh^ 



(B.20) 



The nonhnear Gauss-Seidel smoothing scheme is as follows: 
4sM 



h3 



+ 



2^ 



^+1 



4sM 



(1) -^^ivi / , »+i 



a 1 



li+i I ,,£+1 



> J + ^^i,3 + /i2 



e+1 



1,3 2 



1,3 ' ^« J 



+ 



+ 



(B.21) 



(B.22) 



(B.23) 



At each cell (i, j), we solve for the unknowns Pi^^, ^^U^' '^^ij^ using Cramer's rule. One complete 
smoothing pass is effected when all of the grid points (i, j) have been visited exactly once. The 
smoothing operation is complete when ^max full smoothing passes (or smoothing iterations) have 
been conducted. 
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C Tables 



he 


hf 




4>hf - 4>h, 


2 


Rate 


32 

32 


32 

64 


2.4306 X 10" 


-4 




32 
64 


32 
128 


6.6989 X 10" 


-5 


1.8594 


32 

128 


32 

256 


1.7059 X 10 


-5 


1.9734 


32 
256 


32 
512 


4.2828 X 10" 


-6 


1.9939 



Table 1: Error (successive differences) and calculated convergence rates for the first-order convex 
splitting scheme (3.6) - ( |3.9| ). The parameters are taken to be e = 0.025 (a = 1 — e) and /3 = 0.9, 
and tf = 10. The time step sizes are determined by the quadratic refinement path s = 0.025/i^. 
The global error is predicted to be 0{s) + 0{h?) = 0{h'^), which is confirmed in the test. 



he 


hf 


4>hf - <Phr_ 2 


Rate 


32 
32 


32 
64 


2.5301 X 10""^ 




32 
64 


32 
128 


6.6800 X 10"^ 


1.9213 


32 

128 


32 
256 


1.6818 X 10"^ 


1.9899 


32 
256 


32 
512 


4.2099 X 10"** 


1.9981 



Table 2: Errors (successive differences) and calculated convergence rates for the second-order convex 
splitting scheme (3.12) - ( |3.15| ). The parameters are taken to be e = 0.025 (a = 1 — e), /3 = 0.9, 
and tf = 10. The time steps are determined by the linear refinement path s = 0.05/i. The global 
error is predicted to be 0{s^) + 0{h'^) = 0{h'^), which is confirmed in the test. 
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D Figures 




Multigrid V-Cycle Iterations (m) 



Figure 1: Multigrid convergence test for the second-order convex splitting scheme (3.12) - (3.15). 
The initial conditions are given in (4.1 ); the time step size is held fixed at s = 10 

1, /3 = 0.9, and O = 



while h is allowed 
(0,32)2. 



to vary; and the other parameters are e = 0.025, M = 
^max = 1 in the definition of the smoothing operator (cf. App. B.2), and we conduct exactly 20 
time steps, so that the final time tf = 200. The residuals above are reported only for the last time 
step. We observe that, independent of the space step size h, the norm of the residual is always 
decreasing by about a factor of 5.5 for each multigrid v-cycle iteration. This demonstrates that our 
multigrid scheme is of nearly optimal order. 
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Time 



Figure 2: (Color online.) The computed energies T (3.2, green), T (3.26, red), and F (3.1 blue) 
plotted as functions of time. The computation is effected with our second-order convex-splitting 
scheme ( 3.12[ ) - (3.15), where the initial data are given in (4.1), and M = 1, e = 0.025 (a = 1 — e). 



9. = (0,32) X (0,32), P = 0.01, h 
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and s = 0.1. As pointed out earlier, the energy F is not 



necessarily non-decreasing in time. The only energy that is guaranteed to be non-increasing is 
this property is observed here. But, note that T is nearly identical to J^, which is typically the 



case. 
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Figure 3: The evolution of the density field (j) calculated using the second-order convex-splitting 
scheme and four different time steps: s = 0.01, 1, 10, 20. The parameters are e = 0.025 (q = 1 — e), 
P = 0.9, and h = 1. The plots in the same column have the same pseudo energy J^. 
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Figure 4: The scaled MPFC pseudo energy as a function of time for the evolution from Fig. jsj 
using the fully second order scheme with s = 0.01. Snapshots of the density field 4> are shown in 
Fig. [3| row 1. The numerical and physical parameters for the simulation are given in the caption 
of Fig. |3]and in the text. 
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Figure 5: The evolution of the density field (j) calculated using the first-order convex-splitting scheme 
and four different time steps: s = 0.01, 1, 10, 20. The parameters are e = 0.025, /3 = 0.9 and h = 1. 
The plots in the same column have the same pseudo energy J^. Compare with Fig. [3} 
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Figure 6: Log-log plot of the scaled difference in solutions D{s), defined in (4.2), calculated using 
the second-order convex-splitting scheme. The parameters are tf = 350, /3 = 0.9, e = 0.025 
(a = 1 — e), and h = 1. The baseline calculation is made using the time step size s = si = 0.01; 
and the corresponding baseline evolution of (p is shown in Fig. [3j row 1. 
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Figure 7a: Growth of a polycrystal via freezing of super cooled liquid. The figure shows the density 
field (/) at t = 500. The physical parameters are M = 1, e = 0.25 (a = 1 — e), /3 = 0.9, s = 1 with 
h- 



2048' 
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Figure 7b: Growth of a polycrystal via freezing of super cooled liquid. The figure shows the density 
field (p at t = 1000. The physical parameters are M = 1, e = 0.25 (a = 1 — e), /3 = 0.9, s = 1 with 

2048' 
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Figure 7d: Growth of a polycrystal via freezing of super cooled liquid. The figure shows the density 
field (p at t = 3000. The physical parameters are M = 1, e = 0.25 (a = 1 — e), /3 = 0.9, s = 1 with 
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Figure 8: Growth of a polycrystal via freezing of super cooled liquid. The figure shows the density 
field (f) aX t = 3000. The physical parameters are M = 1, e = 0.25 (q = 1 — e), /3 = 10, s = 1 with 
h = h = . Note the larger value of /3 used here compared with the value used in the simulation 
represented in Figs. [7a| - [7d| Note the near-absence of a grain boundary between the first and 



second grains (about x = 300), compared with the microstructure shown in Fig. 7b 
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Figure 9: (Color online.) Scaled pseudo energy j^j^ as a function of time for the polycrystal 



simulations depicted in Figs. [7a 
/? = 0.9 (Figs. W2 



7d and in Fig. p] The red line corresponds to the parameter 



7d), and the blue line corresponds to /? = 10 (Fig. Here the time axis is on 



the logarithmic scale. 
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Figure 10: An annealed periodic crystal strip, about 20 atom-layers thick, placed in a bath of 
coexistence liquid. For the next test, the results of which are shown in Fig. [TTj we simultaneously 
apply a horizontal shear strain to the bottom row of atoms and pin the top row of atoms. 
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Figure 11: (Color online.) Strain (e^x) fields as functions of the crystal thickness (y-coordinate) at 
times t = 10 (blue), t = 14 (green) and, t = 18 (red). A horizontal shear strain of 4% has been 
applied to the bottom most row of the crystal at t = 0. Simulation results are shown for two cases: 
/3 = 0.2 (triangles, low damping case) and /3 = 9 (squares, high damping case). Note the oscillatory 
nature of the low damping MPFC case (triangles). In this model the crystal rapidly relaxes to the 
strained equilibrium state. The high damping MPFC case behaves much more like the PFC model, 
as expected, with a slow relaxation to equilibrium and no oscillatory overshoot at long times. 
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